//
//  VTKWriter.cpp
//  NSESolver
//
//  Created by Patrick Zulian on 5/9/11.
//  Copyright 2011 ICS - USI - Lugano. All rights reserved.
//


#include "Range.h"
#include "VTKWriter.h"

namespace nsesolver {
    
    VTKWriter::VTKWriter(const std::string & path, const mpi::CartesianTopology &cart)
    : AbstractWriter(path, ".vtk", cart)
    {

    }

    
    void VTKWriter::solverStarted(CFDProblem &problem)
    {
        using std::ofstream;
        using std::string;
        using std::ios;
        using std::ostringstream;
        
        
        const Obstacle & obstacle(problem.getBoundaryConditions().getObstacle());
        Matrix< double > obsMat(obstacle.xmax(), obstacle.ymax());
        obsMat.allSet(0);
        
        Obstacle::ConstObsCellIter it;
        for (it = obstacle.obsCellBegin(); it != obstacle.obsCellEnd(); ++it) {
            obsMat(it->i, it->j) = static_cast< int >(it->type);
        }
        
        writeMatrix(problem, obsMat, "obstacle", "obstacle");
        
        
    }
 
    void VTKWriter::solverStep(CFDProblem &problem, const int q )
    {
        using std::ofstream;
        using std::string;
        using std::ios;
        using std::ostringstream;
        
        
        const Matrixd &p(problem.P()); 
        const Matrixd &u(problem.U()); 
        const Matrixd &v(problem.V());
        
        ostringstream filename;
        filename << path << "/uvp" << q << extension;
        
        ofstream os;
        os.open(filename.str().c_str());
        
        if(os.good()) {
            Matrixd outU(p.xmax() -1, p.ymax() -1);
            Matrixd outV(p.xmax() -1, p.ymax() -1);
            Matrixd outP(p.xmax() -1, p.ymax() -1);

            AbstractWriter::auxShift(outU, u, 1, 0);
            AbstractWriter::auxShift(outV, v, 0, 1);
            AbstractWriter::auxShift(outP, p, 1, 1);
            
            
            Range  rangeX(0, outU.xmax() - !isBoundary(East, cart) ), 
                   rangeY(0, outU.ymax() - !isBoundary(North, cart) );
            
            
            
            const MatrixHandle< Matrixd > & uh(outU(rangeX, rangeY)), 
            vh(outV(rangeX, rangeY)), 
            ph(outP(rangeX, rangeY));
            
            ostringstream vtkHeader;
            
            vtkHeader << "# vtk DataFile Version 2.0\n";
            vtkHeader << "Velocities and Pressure NSE\n";
            vtkHeader << "ASCII\n\n";
            vtkHeader << "DATASET STRUCTURED_POINTS\n";
            vtkHeader << "DIMENSIONS " << rangeX.getExtent() << " " << rangeY.getExtent() << " 1\n";
            vtkHeader << "ORIGIN " << problem.IGlobal() << " " << problem.JGlobal() << " 0.0\n";
            vtkHeader << "SPACING " << problem.DX() << " " << problem.DY() << " 0.0\n\n";            
           
            vtkHeader << "POINT_DATA " << rangeX.getExtent() * rangeY.getExtent() << "\n";
            vtkHeader << "SCALARS pressure double\n";
            vtkHeader << "LOOKUP_TABLE default\n";
            os << vtkHeader.str();
            
            os << ph;
            
            
            os << "\nVECTORS velocity double\n";
            
            for (int y = 0; y < uh.ymax(); ++y) {
                for (int x = 0; x < uh.xmax(); ++x) {
                    os << uh(x, y) << " " << vh(x, y) << " 0.0\n";
                }
            }
        }
        os.close();
       
    }
    
    
    void VTKWriter::writeMatrix(const CFDProblem &problem, const Matrixd &m, const std::string &varName, const std::string & filename ) const
    {
        using std::ofstream;
        using std::string;
        using std::ios;
        using std::ostringstream;
        
    
        ostringstream filePath;
        filePath << path << "/" << filename << extension;
        
        ofstream os;
        os.open(filePath.str().c_str());
        
        if(os.good()) {
            
            ostringstream vtkHeader;
            
            vtkHeader << "# vtk DataFile Version 2.0\n";
            vtkHeader << varName << " of NSE\n";
            vtkHeader << "ASCII\n\n";
            vtkHeader << "DATASET STRUCTURED_POINTS\n";
            vtkHeader << "DIMENSIONS " << m.xmax() << " " << m.ymax() << " 1\n";
            vtkHeader << "ORIGIN " << problem.IGlobal() << " " << problem.JGlobal() << " 0.0\n";
            vtkHeader << "SPACING " << problem.DX() << " " << problem.DY() << " 0.0\n\n";            
            
            vtkHeader << "POINT_DATA " << m.size() << "\n";
            vtkHeader << "SCALARS " << varName << " double\n";
            vtkHeader << "LOOKUP_TABLE default\n";
            os << vtkHeader.str();
            
            os << m;
            
            
            
        }
        os.close();
    }
    

    
}
